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ABSTRACT 

Planets migrate due to the recoil they experience from scattering solid (planetesimal) bodies. To first order, 
the torques exerted by the interior and exterior disks cancel, analogous to the cancellation of the torques from 
the gravitational interaction with the gas (type I migration). Assuming the dispersion-dominated regime and 
power-laws characterized by indices a and j3 for the surface density and eccentricity profiles, we calculate the 
net torque on the planet. We consider both distant encounters and close (orbit-crossing) encounters. We find 
that the close and distant encounter torques have opposite signs with respect to their a and ft dependences; 
and that the torque is especially sensitive to the eccentricity gradient {ft). Compared to type-I migration due to 
excitation of density waves, the planetesimal-driven migration rate is generally lower due to the lower surface 
density of solids in gas-rich disk, although this may be partially or fully offset when their eccentricity and 
inclination are small. Allowing for the feedback of the planet on the planetesimal disk through viscous stirring, 
we find that under certain conditions a self-regulated migration scenario emerges, in which the planet migrates 
at a steady pace that approaches the rate corresponding to the one-sided torque. If the local planetesimal disk 
mass to planet mass ratio is low, however, migration stalls. We quantify the boundaries separating the three 
migration regimes. 



Keywords: planetary systems: protoplanetary disks - 
ods: analytical 



1. INTRODUCTION 

Bodies immersed in gaseous or particle disks migrate ra- 
dially. Very small particles, strongly coupled to the gas, are 
carried by the gas. Thus, they fol low the accretion fl ow or are 
dispersed by turbulent motions (Ciesla 20091 120101) . Larger 
particles tend to move on Keplerian orbits. However, in proto- 
planetary disks the gas is partially pressure-supported, which 
causes soli ds to drift inwards due to the headwind th ey ex- 
periences (Ada chi et al J [19761 IWeTd enschillirigl [T977h . This 
effect peaks for ~m-size bodies (or their aerodynamic equiv- 
alents) at which they spiral in in as little as ~100 orbital pe- 
riods. Larger, km-size bodies (planetesimals) are more re- 
sistant against drag-induced orbital decay due to their large 
inertia. The motions of these bodies will be predominantly 
determined by gravitational encounters, rather than gas drag. 

The gravitational interaction with the gas also causes a drag 
force on the planet. The picture here is that of a massive body 
gravitationally perturbing the disks, which causes an excess 
density structure that backreacts on the planet. One can re- 
gard the force that the planets experiences a manifestation of 
dynamical friction - a concept that is perhaps more famil- 
iar with collisionle ss systems, but which can also be applied 
to gaseous disks (lOstrikeri [l999t iKim & Kiml 120071 l2009t 
iMutoet alJl20lU iLee & Stahlerll2()l 11) . When the planet is 
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small, the resulting migration from th e gravitational interac- 
tion with the gas i s known as type I dGoldreich & Tremaine 
Il980t IWard1ll986l) . Like dynamical friction, the type-I mi- 
gration rate increases linearly with mass. Although rather 
insignificant for planetesimals, it becomes very efficient for 
Earth-mass planets resulting in migration timescales as short 
as 10 s yr at 1 AU dTanaka et al.ll2002l) . 

For these reasons (gas-driven) migration is often invoked 
to explain the existence of close-in, Neptune- and Jupiter- 
mass planets ('hot Jupiters'), since conditions very close to 
the star are thought to b e ill-suited to form giant planets in 
situ dlkoma &H ori 2012). However, a clear understanding of 
type-I migration is somewhat complicated by the fact that it 
is a higher order effect; that is, the net torque on the planet 
results from a near cancellation of two large but opposite 
torques, corresponding to the respective contributions from 
the inner and outer disks. In addition, the net co-orbital and 
Lindblad torque may have different signs. As a result the sign 
of type-I migration is very sensitive to the local distribution 
of matter, which in turn is determined by the thermodynam ic 
properties of the disk (e.g., Paardekooper & Mellema 2006). 

Similar to 'gas-driven' migration, scattering of solid bodies 
also causes a planet to migrate. This effect of planetesimal- 
driven migration (PPM) has been mostly explored through 
Af-body studies (lHahn &Malhotral 1 1 999t iKirshet all 120091: 
iBromley & Kenyonll201 lbtlCapobianco et al.ll2.01 ll) . In some 
cases, these authors found an migration instability, at which 
the planet migrates at a rate determined by the one-sided 
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torque (llda et al.ll2000l) . Under these conditions, PDM is fast. 

Other studies have investi gated the embryo-planetesimal 
interaction analytically (e.g. Jldalll990t llda & Makinolll993t 
iTanaka & Idalll996t lRafikovN2003al) Mostly, these studies 
consider the effect of the embryo on the planetesimal disk, 
e.g., the rate at which the protoplanet excites the planetesi- 
mal's eccentricity or how it opens a gap by scattering. 

In this paper, on the other hand, we will study the recoil of 
the scattering on the planet for given planetesimal properties. 
These calculations provide, for the first time, an analytic ex- 
pression for the two-sided torque for planetesimal scattering 
- the analogue to the type-I migration torque. 

We assume the following: (i) a smooth disk where the spa- 
tial distribution of surface density and eccentricity are power- 
laws; (ii) the dispersion-dominated regime (relative velocities 
are given by the eccentricity of the planetesimals at close en- 
counter); (iii) Keplerian orbits for the planetesimals; (iv) a 
circular orbit for the planet. We account for both distant and 
close encounters, corresponding to orbits that do or do not 
cross the planet (see Fig . [TJ . We then compute the recoil of 
the planet due to scatterings with planetesimals on orbits both 
interior and exterior to the planet, which results in the PDM 
rate. 

PDM can be divided into three regimes, depending on the 
ratio of the planet mass compared to the mass of the solids 
with which it interacts: 

1. Low mass planets. They do not exert a (strong) feed- 
back on the disks. Correspondingly, the gradients 
in planetesimal's eccentricity and surface density are 
those of the background disk {a and /? in Fig.[TJ and 
can be assumed fixed during the migration; 

2. Massive planets. They have difficulty to migrate over 
large distances due to their inertia. Instead, t he planet 
scatter s away the planetesimals, leaving a gap (Rafikov 

l2003allbh . 

3. Intermediate-mass planets. They exert some feedback 
on the disk but not enough to halt their migration. 

In § |2f{3] the first regime is assumed. In §|2] the calculation 
for the migration rate due to distant and close encounters are 
presented. In § [3] the net torque and the corresponding mi- 
gration timescale are computed and compared to the type-I 
migration timescale. Furthermore, the approach is sketched 
how a distribution in eccentricity must be incorporated. In §|4] 
the importance of diffusive motions ('noise') is investigated. 
Then, in § we focus on the third regime and find that the 
migrating planet regulates the local eccentricity profile. Fur- 
thermore, we will outline the boundaries dividing the regimes 
and find that the intermediate regime covers a large region of 
the parameter space. We summarize our results in §|6] 

2. CALCULATION OF THE MIGRATION RATE 

2.1. Statement of the problem and methodology 

We consider the following setup. A planet of mass M p 
moves on a circular, non-inclined orbit at a reference disk 
radius oo m the equatorial plane. The planet interacts grav- 
itationally with planetesimals of mass m <K M p that are char- 
acterized by standard Keplerian orbital elements: semi-major 
axis a, inclination i, eccentricity e, and phase angles. The 
surface density of the planetesimals is given by X(a), and the 
planetesimals are assumed to be randomly distributed in their 
phase angles: mean anomaly f, and argument of periapsis, a>. 



Distant 
(Repulsive) 


1 

i Close 
i | (Attractive) 
i 


1 

i Distant 

| (Repulsive)^^ 

i ^^^^ 




i 

1 ^^"^ 

\ j t e o"o 


i 
i 




: ' ^in^in ^ou^oli 


1 




. ),( J t 

J_|fin___Jfl 

1 

: 1 


*l p 

E 0U 

I 
1 

1 




aim ao 
Disk semi major axis, a 





Figure 1. Sketch of the disk profile. A planet on a circular orbit (e = 0) 
at semi-major axis ag interacts with planetesimals either through close or 
distant encounters. Planetesimals that are able to cross the planet's semi- 
major axis interact via close encounters; otherwise the encounters are distant. 
We allow for a power-lawprofile of surface density and eccentricity with 
indices a, jj (see Equation (T)) and compute the net migration rate dag /dt that 
the planet experiences due to scattering of the planetesimals in the dispersion- 
dominated regime. A nonzero p causes the transition between the regimes 
(dotted and dashed vertical lines) to shift by an amount xfieiaa (see text). 

The calculations allow for gradients in 2 and e; specifically, 
they are assumed to be a power-laws with indices a and /3: 

S(a) = S (-) ; «(a) = «o(-) , (D 
\ao/ \ao) 

where So and eo are reference values that correspond to the 
surface density and eccentricity at the semi-major axis of the 
planet (a = a®). 

This configuration is sketched in Fig.Q] where the surface 
density and eccentricity follow a power-law profile. There are 
two ways in which the planetesimals interact with the planet 
- close and distant encounters - dependent on whether or not 
they are able to cross the planet's orbit. In close encounters, 
the planet tends to scatter planetesimals from the exterior disk 
to the interior disk and vice versa. This is a dispersive process: 
the net separati on after the scattering on average increases 
(Rafikov 2003a). But due to the recoil from the scattering, 
the planet moves in the direction of the denser planetesimal 
belt. Thus, from its perspective it is attracted by the belt, al- 
though it may strongly excite the belt over the course of its 
migration. 

Distant encounters are repulsive, in the sense that they push 
the planet away from a planetesimal belt. Since the planet is 
assumed to move on a circular orbit, the boundary between 
the distant and close encounter regimes is determined by the 
periapsis (a(l - e[a])) and apoapsis (a(l + e[a])) of the plan- 
etesimals' elliptical orbits. Planetesimals of semi-major axes 
less then a; n or larger than a ou interact through distant encoun- 
ters; planetesimals of semi-axes a m < a < a ou through close 
encounters (see Fig.[TJ. Here, a m , a ou are given by: 

Oinlou [1 + efainlou)] = «o, (2) 

where the upper sign corresponds to a ou , the lower to a m . 
When Equation ([1) is linearlized by writing e as eo + 
/teo(flin|ou _ fl o)/ao> we can so l ve for ^iniou to obtain 

flin|ou ~ a ± e a + e^fia^. (3) 

A positive /? (depicted in Fig.[T|) therefore shifts the transition 
between the distant and the close encounter region to larger a 
by an amount ?*f}e\a§. 
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Table 1 

List of frequently-used symbols. 



Symbol Description 



Avj Change in the i-th component of the relative velocity 

T Dimensional torque on planet 

A Coulomb factor 

X(a), Zo Surface density of planetesimals at disk radius a or at reference 
radius ao 

Clo, fl a Orbital frequency at corresponding to the reference position or 

to semi-major axis a 

a Exponent in surface density power-law (Equation (TJ) 

ji Exponent in eccentricity power-law (Equation (TJ) 

yx Dimension less torque for close or distant encounters 

(Equation J24l ) 

y cv Curvature component of dimensionless torque (Equation J24J) 

7v Gradient component of dimensionless torque (Equation i24\ ) 

S Exponent in inclination power-law 

(p Azimuthal coordinate in cylindrical coordinate system 

v Viscosity or diffusion rate in semi-major axis (Equation {43)) 

8 True anomaly (0 = indicates periapsis; Appendix) 

D[Avj] Diffusion coe fficient: rate of change in relative velocity 

(Equation j[6i ) 

Gjv Newton's gravitational constant 

M+ Stellar mass 

M p Planet mass 

Qpi Combination of qj and g p , Equation 457t 

R Radial coordinate in cylindrical units 

R], Hill radius Equation (4) 

Pr z Probability dens ity o f finding a planetesimal near R = ao and 

z = (Equation flSH ) 

Tmigr Timescale to migrate globally over distance ao 

^inigr Timescale to migrate locally over distance egao 

T syn Synodical period 

^typc-i Type-I migration timescale 

T vs Viscous stirring timescale (Equation \49\ ) 

Vko Kepler (orbital) velocity corresponding to ao 

"[in.ou] Inner or outer-most semi-major axis from where planetesimals 

cross the planet's orbit 

ao Semi-major axis of the planet; reference radius 

b Distance between semimajor axis planet and planetesimal 

[e, eo\ Eccentricity of planetesimals (at a = ao) 

ei, Hill eccentricity (Equation {5)) 

e* Lower range of the Hill ecc entri city for which self-regulated 

migration applies (Equation )59E ) 

(' Inclination of planetesimals 

/a Coulomb term, f\ = log(l + A 2 ) 

g\ Coulomb term, g^ = A 2 /(l + A 2 ) 

m Mass of individual planetesimal 

gj Dimensionless 'disk mass' (Equation )10> ) 

q„ Dimensionless mass of the planet (=M p /M+ ) 

r Radial coordinate in polar coordinate system 

v Relative velocity between planet and planetesimal at a = ao 

Vj Relative velocity of f'-th component 

z Vertical coordinate in cylindrical units 



For (very) low e the distinction between close and distant 
encounters is no longer determined by the eccentricity but by 
the Keplerian shear of the disk. In the limit of zero eccentric- 
ity and inclination, the b order between distant and close en- 
counters lies at ~2.5R h (lNishidalll983l: iPetit & Henonll 19861) 
where R/ t is the Hill radius, 



Rh 



ao 



M p 
3M* 



1/3 



«() 



(4) 



with M p the planet's mass, the stellar mass, and q p = 
Mp/Mi,. Furthermore, in the zero-eccentricity limit planetes- 
imals of semi-major axis similar to the planet travel on horse- 
shoe orbits and do not enter the Hill sphere. Interactions in 
this shear-dominated regime are qualitatively different from 
the high velocity regime, where random motions (caused by 
the eccentricity of the planetesimal) dominate. Random mo- 



tions start to dominate over shear motions when eVx k. RhQ) 
or for e > (M p /3M i ,) 1/i , where V^o = ao^h is the Keplerian 
orbital velocity and ao and Qo the local orbital frequency. It 
is sometimes convenient to express eccentricities in terms of 
RhClo rather than V^; i.e., 



so called Hill eccentricities. In this work it is assumed that the 
dispersion dominated regime applies: e% > 1. 

2.2. The contribution f mm distant encounters 

lHasegawa & N akazawa d 19901) have calculated the change 
in relative semi-maj or axis due to a distan t encounter among 
two bodies (see also Hen on & PetitHl986l) . Assuming a uni- 
form distribution of phase angles, the average change for the 
encounter becomes: 



<Afc> = 



b 5 



where 



C 



54 [2tfo(2/3) + *i(2/3)] 



30.1 



(6) 



(7) 



is a constant, b - a - the separation in semi-major axis be- 
tween planet and planetesimal, and K v is the modified Bessel 
function of the second kind of order v. The encounter is 
always repulsive; (Ab) has the same sign as b and, after 
phase-averaging, is independent of eccentricity and inclina- 
tion. When we consider the interacting bodies to be a plan- 
etesimal of mass m and a planet of mass M p » m, the plan- 
etesimal will experience the largest change in its semi-major 
axis. But due to the recoil effect, the planet experiences a 
change of 



(Ab) M = - 



m + M, 



■(Ab) * -—(Ab). 

Mr, 



(8) 



The rate at which a planet migrates due to encounters with 
planetesimals at distance b is given by Equation © times the 
encounter rate. To first order the encounter rate for an impact 
parameter b is ||/7|Q ( )Eo/w, where we took the local value of 
the surface density at the planet's position and linearized the 
(Keplerian) shear. If the planet interacts only with planetesi- 
mals on one side of its orbits, e.g., with planetesimals exterior 
to it (b > 0), the migration rate becomes: 

da\ 1 r°° 3 

— =-tt db-Z b£l (Ab) 

CYMpalQ.0 



1.67-^(00^0), 



(9) 



(when the inner disk is considered, the sign will be positive) 
where the subscript 'ls-di' refers to 'one-sided' and 'distant 
encounters'. In Equation (0 the dimensionless disk mass qj 
is defined as: 



qd(a) 



S(a)fl 2 
M+ 



10" 



2, 



10 g cm 



\2+a 



a i _ 



(10) 

where a\ is a reference radius (say 1 AU) and Ei the surface 
density at a = a\. Note that is a local quantity that depends 
on disk radius ao. In the outer disk, qj may be significantly 
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larger, since ices will contribute to the solid fraction and 2 + a 
is typically positive. 

Equation (O gives the migration rate due to distant encoun- 
ters for one side of the disk. This zeroth order effect scales as 
oce Q 3 . As the contribution from the interior disk will have the 
opposite sign, these terms cancel to zeroth order. 

Therefore, we consider the next order contributions. These 
arise due to gradients in eccentricity and surface density 
(Fig- Q3 and due to higher order approximation of the velocity 
field around the planet instead of the sheering sheet approx- 
imation used in the local Hill formalism which Equation (O 
relies on. These latter effects are referred to as curvature ef- 
fects. To obtain the higher order term , we have used the for- 
malism of linear density wave theory dGoldreich & Tremaind 
ll980tlWardiri986lll997h . These provide us with a formula for 
the torque density - the torque per unit disk radius - the planet 
exerts on the disk. In AppendixlAlwe derive the expressions 
for the torque that the planet experiences. The result is: 



Tis- 



fO.84 



-qdq P - 



(11) 



M p (a<A)) 2 

r 2s _ di -5.7 + 2.5(-a + 2(3) 

MJa Q. Y ei 



In these expressions Ti s denotes the torque the planet experi- 
ences due to one side of the disk only, where the upper sign 
corresponds to the exterior disk (integration over the distant 
encounter region where b is positive) and the lower sign to 
the interior disk. The two-sided torque T2 S corresponds to 
contributions from both sides of the disk. As remarked, this 
expression is an order higher in eo than the one-sided torque 
- but still significant. 

The torque adds or removes angular momentum to the 
planet at a rate dljdt where Z- = o^M p Qq. Since dZ/df = 
M p (dfl/df)d(flgQo)/dfl we obtain the migration rate as 



2Y 



dao 



and it can be verified that Equation ( TTTb is consistent with 
Equation When accounting for both sides of the disks, 
the migration rate becomes 



dao\ 



n.3 + 5.0(-a + 2/3)q d q P/ 
5 (fl(A 



i). (14) 



The sign of the migration thus depends on the values of a 
and (3. A large, positive value of a implies that interactions 
with the exterior disk will dominate, which pushes the planet 
inwards. Positive /3 implies that a m lies closer (in absolute 
terms) to ao than a ou , which tilts the balance in favor of the 
inner disk. However, due to the large negative value of the 
curvature term, the direction of migration tends to be inwards 
in most cases. Finally, lower eo increases the importance of 
distant encounters as both a m and a ou move closer to ao. 

Distant encounters represent only one side of the medal. 
Close encounters reverse the sign and have the opposite de- 
pendences on a and p. For the net migration rate both must 
be considered. 

2.3. The contribution from close encounters 

A scattering of 2 bodies rotates the relative velocity vector 
v, while preserving its absolute value. For the migration rate 
it is the change in the azimuthal component of v, Au^, that 



matters and this component is generally not conserved. To- 
gether with the encounter rate they determine the force that 
th e planet experiences . 

iBinney & Tremainel (120081) have calculated the diffusion 
coefficients - the rate of change in the components of v. Af- 
ter integration over impact paramet ers they obtain the rate of 
change in the parallel component ( Binn ev & Tremainell2008l 
their Eqs. L.ll): 



D[Au||] = 2nnv 



G 2 N m(M p + rri) 



/a- 



(15) 



where /a = ln(l + A 2 ) is a Coulomb term which is assumed 
to be constant (i.e., independent of velocity) in the following, 
Gn Newton's gravitational constant, m the mass of the field 
particles (planetesimals), n the local number density (at the 
planet's position), and v = |v] the relative velocity between 
the planet and the unperturbed (Keplerian) orbit of the plan- 
etesim als at the interaction point (sometimes called collision 
orbits; Tanaka & Idal ll996l) . 

Equation ( [TBI ) does not include the effects of the solar grav- 
ity, which can be effective t o change the o r bital o f planetes- 
imals during the scattering. Tana ka & Idal ([1996) examined 
the effects of the solar gravity and found that the effect of the 
solar gravity cancels, after averaging over the (uniformly dis- 
tributed) phase angles (i.e., the longitudes of periapsis and as- 
cending node). This cancellation is related to the fact that the 
unperturbed (i.e., Keplerian) axisymmetric particulate disk 
does not exert any torque on the planet. Thus the effect of 
the solar gravity during scattering will cancel even in our case 
where non-local, i.e., curvature terms, effects are included. 

For the individua l components we have 
dBinnev & Tremaindl2008l their Equation [7.89]): 



D[Avj] = jD[Av n ]. 



(16) 



(jj) Here we consider the change in the azimuthal (0) component 



and further assume that M p » m: 



D[Av 4 ,]=2nGlnM p mf h ^-. 



(17) 



Equation ( TTTb is nothing else than the azimuthal force exerted 
on the protoplanet due to dynamical friction with the planetes- 
imals. 

However, Equation ( fT71 ) is only valid when the density n 
and velocity field v are uniform. This is certainly not the case 
in a Keplerian disk; particles of different semi-major axis a 
will have a different (relative) velocity at the point where they 
interact with the planet, that is, at a — ao. Equation (TTTb has 
to be convolved over the semi-major axis. The same holds 
for the density, n. For particles traveling on a Kepler orbit 
with eccentricity e and semi-major axis a, the density is not 
constant as the particle's velocity depends on its position (true 
anomaly). 

Let us introduce the projection operator Ps^, defined such 
that Pr zt i,dR(Rd<l>)dz gives the probability of finding a particle 
with orbital elements a, e, i and random phase angles in the 
interval [R, R + dR; (f> + d</>; z + dz], where (R, <p, z) are cylindri- 
cal coordinates. Clearly, P Rz</l is a function of the properties of 
the particle (a, e, i) as well as the position at which it is eval- 
uated, as given by the coordinates (R,z, <f>). For the latter we 
assume azimuthal symmetry, R - ao, and z — corresponding 
to the position of the planet and define a new, more specific, 
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projection operator: 

pR Z (a,i,e) — 2naoPR < p z (a,i,e;R = ao,(p,z = 0), (18) 

such that pR : dRdz gives the probability that an (a, e, /)-particle 
can be found in the equatorial plane at ao at arbitrary (p. 

Using Equation ( fT8l we obtain the contribution to the den- 
sity n from particles of semi-major axis a. Since the total mass 
of particles in [a, a + da] equals dm - 2ndL(a)da: 

1 1 ~ aL{a) 

dn — — 2jTaY.(a)PR <p ,(a; ao, z — 0)da — — Pr, da. 

m m ~ ao 

(19) 

With this notation the specific force on the protoplanet, 
Equation ( fT71 ), becomes: 



F = D[Atv] = 2nG 2 N M p f A 



f 

Ja\ n 



daP Rz Z(a)\- \% (20) 
ao v i 



where all quantities in the integrand are functions of semi- 
major axis, a. The integration proceeds over the close en- 
counter region. 

The vertical component of the torque exerted on the proto- 
planet is given by T z = aoM p Ff 



T c i = F^aoM p = 2nG N M p aofK 



-Ms ? 



da P fe Eo — 

\ao) y J 

(21) 

where we have inserted Equation (0 for 2(a). This integral 
gives the migration rate due to close encounters. To solve 
it, the velocity field of the planetesimals near the protoplanet 
(the v and v$ terms) and Pr- must be expressed as function 
of a. These steps are outlined in AppendixlBl Equation (fJTJ 
also depends on the inclination of the particles - a thinner disk 
will, for example, increases the density of particles (so that P 
increases) and additionally decreases the relative velocity v 
(since v z decreases). These effects increase the magnitude of 
T c i and therefore the migration rate. 

In AppendixlBl we perform the calculations for the case 
where i = e/2 — the equilibrium solution- and a case where 
i <s e. We obtain, for the one-sided torques: 



(a £lo) 2 M p 



Aqdq P 



±1.1 
+ 1.3 



(i = e/2 «: 1); 



(i «: e «: 1); 



(22) 



(the upper sign corresponds to the torque that the exterior disk 
exerts on the planet; the lower to the interior disk); and for the 
two-sided torque: 



2s-cl 



(a<A)) 2 M„ 



f\qdq P 



-0.7 + 2.0(a - 3/3) 



1.6 + 2.3(a-2/3-5) 



(i = e/2); 



(i <k e); 



(23) 

where 6 is the exponent of the inclination dependence with 
disk radius, i.e., the inclination equivalent of [3. Contrary to 
the distant encounter torque (Equation Sl2\ ). T2 S - C \ increases 
with increasing a, since for close encounter the planet tends 
to move in the direction of the more massive planetesimal 
belt. In addition, Equation ( f23l displays a negative depen- 
dence on /3, which can be understood since the cross section 
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Figure 2. Sign of migration due to scattering by close and distant encounters 
for the two-sided torque with (' = e/2. We plot the line in the (a,fl) plane (re- 
ferring to the indices in surface density and eccentricity) where the migration 
direction changes sign. For close encounters migration is outwards below 
the blue, dashed line. When distant encounters are included (black lines) the 
dividing line shifts down, dependent on the value of the Coulomb factor /a ■ 

for encounters is largest when the relative velocity (eccentric- 
ity) is lowest. More generally, the relative importances of the 
gradient terms follow directly from Equation (|2TT> (or even 
Equation (15[). In it Pr z reflects the scaleheight dependence 
and consequently contributes a term -6; S(a) a term a; and 
ity/o 3 , which is proportional to e~ 2 , a term -2/3. Note the dif- 
ference in the curvature term (the first term in Equation ( f23l ) 
between the equilibrium case (<? = i/2) and the thin disk case 
(i <K <?): in the latter it is positive, whereas in the former it is 
negative. However, the curvature term for distant encounters 
(Equation Sl2\ ) is always negative. 

Thus, the planet is attracted towards quiescent and 
mas sive planetesimal belts, in line with previous stud- 
ies (iTanaka & Idalll999t IPavne et alj|2009t ICapobianco et al l 
1201 ll) . During its migration, the planet scatters away many 
planetesimals, which in turn affects their spatial and dy- 
namical distribution. The amount with which the scattering 
changes the (effective) values of a and /3, depends on the rel- 
ative masses of the planetesimal belt and the planet (see §|5]l. 

In the remainder, we will focus on the equilibrium solution 
(i/e = 2 and 5-/3) since thi s is the expect ed ratio for the 
dispersion-dominated regime (llda et al.lll993l) . We consider 
an extension to an eccentricity distribution in § 13.41 

3. TORQUES AND MIGRATION RATES 

3.1. Outwards or inwards 

Let us decompose the dimensional torque Tx for an interac- 
tion X as follows: 



M p (a Oo)- 



■ = qdq p F x (eo, k) x Yx(a,P, 6) 



(24) 



= qdq P F x (eo, io) x [% v + fy(a + g(J3, 6))] (25) 

Here, Fx(eo, io) is a function of eo and io only (without a 
numerical prefactor) and jx is the dimensionless torque for 
interaction X. For a two sided toques, y is further decom- 
posed into a curvature term % v and a gradient term yy, de- 
fined such that it is proportional to a in yx- For example, 
for the two-sided, close encounter torque of Equation d23l ), 

r^-cf = ("°- 7 + 2 -°[ ff - Wa> F = e Q 2 , % v = -0.7/,, 
yv = 2.0/, and a = -3/3. We have compiled a list of dimen- 
sionless torques in Table [2] 
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Table 2 

Dimensionless torques for planetesimal-driven migration y. 



Torque type 
(1) 




Symbol 

yx 

(2) 


Leading term 

F(eo,k) 
(3) 


Zeroth-order term 
(4) 


Curvature term 
(5) 


rv 

(6) 


Gradient terms 
(V) 


Distant, 1 -sided 
Distant, 2-sided 
Close, 1 -sided, i 
Close, 1 -sided, i 
Close, 2-sided, i 
Close, 2-sided, i 


= e/2 
<K e 
= e/2 
<K e 


yis-cu 

r2s-di 
Ji=el%) 
' ls-cl 

(/<«') 

''ls-cl 
J?=efa 
' 2s-cl 

^2s-cl 


% 3 

r [ e - 2 

e- 2 

■-1 -2 

'o e o 


+0.836 

±1.14/ A 
±1.27/ A 


-5.66 

-0.66/a 
1.63/a 


-2.51 

1.97/a 
2.28/a 


(a - 2(3) 

(a - 3/3) 
(a -2/3- S) 



Note. — Summary of dimensionless torques expression for planetesimal scattering in the dispersion-dominated regime (see Equation i24i). 
Column (1): torque: distant or close; one or two sided; and the inclination model. Column (2): corresponding symbol. Column (3): the order of 
the contribution to the torque in terms of the eccentricity and inclination at the reference position (see Equation J24b). Column (4): the zeroth 
order contribution from one side of the disk (the upper sign corresponds to the outer disk; the lower to the inner); Column (5): the contribution 
to y that arises due to curvature, i.e., the deviation from the linear approximation; Column (6) and (7): the contribution to y due to gradients in 
surface density and eccentricity (see Fig.[TJ. 



The direction of the migration is determined by the sign of 
Ttot = Tci + Tdi and is positive for outwards migration, negative 
for inwards migration. It depends on a, /3 and on the value of 
the Coulomb term f A = log(l + A 2 ). For example: 



-6.3 - 0.5a - 0.9/3 (/ A = 1) 
-7.6 + 3.4a- 12.7/3 (/a = 3) 
-8.9 + 7.3a -24/? (f A = 5) 



(26) 



Thus, for increasing f A , |y tot ] generally becomes larger with 
the sign of the migration more likely to be determined by 
that of the close encounter contribution. This is illustrated 
in Fig. |2] For the range in a and /3 displayed in Fig. [2] y& 
is always negative (inward migration), mainly due to the large 
(negative) value of the curvature term. Close encounters more 
readily give rise to outward migration, but require a more 
massive outer disk (large a) and/or a sufficiently low eccen- 
tricity gradient (reflecting a dynamically colder state with 
which the planet interacts more strongly). Accounting for 
both distant and close interactions shifts the boundary line di- 
viding the inward and outward migration regime down. The 
amount of the shift depends on the Coulomb parameter f A , 
which therefore translates in a relative measure of the impor- 
tance of close vs. distant interactions. 

What is the expected value of /a = log(l + A 2 )? Here, 
A - b max /bgo is the ratio for the l argest and typical im- 
pact radii (Bin nev & Tremainel 120081) . For the former we 
may substitute the disk scaleheight, aoi'o while the latter - 
the impact radius that causes a ji/2 change in the relative ve- 
locity after the scattering - is, in the high velocity regime, 
bgo = GNM p /(eVko) 2 = 3Rh/e\ where e/, is the Hill eccen- 
tricity (Equation (0). Assuming the dispersion-dominated 
regime, i = e/2, f A log(l + e^/6). For e h > 3, f A > 3 
and the close encounter contribution typically determines the 
sign. But at small Hill eccentricity f A - 1 and distant encoun- 
ters become more important (see Fig.|2]i. 

Assuming that the random motion of planetesimals is bal- 
anced by gas drag, it follows that the Hill eccentricity e;, is 
independent of the planet mass, and has a weak (ozX 1 ^ 5 ) de- 
pendence on the planetesimal radius and gas density (smaller 
planetesimals or denser gas results in lower e/,), disk radius 
(larger qp have lower eh). A typical range may be e/, 3-8 
(Kokubo & Ida 2002). When the gas is absent, e/, will in- 
crease with time until it is equilibrated by collisional damp- 
ing. 



3.2. The migration timescale 
The migration timescale is defined: 



1 dap 

1 micr — I j 



2Tv 



v M ; ,fl 2 Q ( )J 



(27) 



In terms of y tot (Equation (f24|) the timescale corresponding 
to the two-sided, i — e/2 torque becomes: 



-On 



2ytot&#p~ 

i— ) 2 

V0.02/ 



(28) 



:2X 10 3 



ytot 



10 



lio- 4 J lio- 6 J ■ 



Alternatively, the migration timescale can be expressed in 
terms of Hill eccentricity e/, using Equation (0: 



T = 

1 misr 



3 2/3 ytot<7</<7 P 

(f) 



1/3 ^0 

2 



-I 



(29) 



;2.2x W 



ytot 



10 



(io^) (io^) °° 



which is useful since the expressions are valid only for e/, > 1 . 
For reference, we also give the one-sided migration timescale 
due to close encounters corresponding to yfcf.j : 



1 



2.2/a q d q p 
l.lxl0 4 (=| 



(30) 



(0.02) do" 4 ) do" 6 ) 



3.3. Comparison with type I migration 

The migra t ion tim escale for type I migration is given by 
lTanakaetalJd2002l) : 



1 



type-I ■ 



2y/0) \ V k o 



;,o«o 



Mp 



(31) 



2 E quation {30} is consistent with Equat ions (19) and (2 0) of llda et al] 
120001) . On the other hand, Equation (23) of llda et all 12000) is inconsistent 
with Equation 430> by a large factor (~10-100), since several numerical con- 
stants were omitted. 
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-0.2 - 



ft(ff«o) 

- - P„(a= Uovo) 



0.0 



0.5 



1.0 



1.5 2.0 

ef oio 



2.5 



3.0 



3.5 



Figure 3. Eccentricity distributions: (black curve) Rayleigh distribution for 
<t = a-go; (dashed curve) Rayleigh distribution for tr = \.\<r e o, which corre- 
sponds to the distribution at separation b = 0.1/jS in the case where cr e is a 
power-law with exponent ft; (red curve) correction distribution pbPc ; (gray 
curve) approximation to Pn(cT[b]) in terms of P/((<x f o) and Pc(°"eo)- 



:2.5 x 10 4 



yi(a)\ 1 j c s /V k 
0.05 



\ 2 / 1g \ X ( 1p y 

J uo- 2 i \ io 6 ^ 



■1 



where y/ denotes the dimensionless torque for type-I migra- 
tion, and q g = Y lg a 2 ) /M e the dimensionless disk mass in gas. 
The expression of ji depends on several factors, including the 
type of torque considered (co-orbital, Lindblad), the surface 
density exponent, the pr essure expon e nt, an d the exponent 
for the scale height ("see | Tanaka et ail (12002|) and successor 
works, e.g. , [D' Angelo & Lubowl 120 lOt iPaa rdekoop er et al.l 
120101 1201 11 iMassefcoi lb . Generally, the situation for type-I 
is more complex, because of the pressure effects and the en- 
ergy transfer in gaseous disks. Nevertheless, the similarity 
between Equations ( l28l l and (IBTT l is striking. In fact Equa- 
tions ( l28l l and ( |3TT > are the same, except that the eccentricity 
in Equation OTb has been substituted by cJViX) an d the sur- 
face density in solids (So) by that of the gas (E^o)- 

3.4. The effect of an eccentricity distribution 

Up till now we have neglected an intrinsic probability distri- 
bution in eccentricity and inclination, w hich may be expected 
from gravitationally-int eracting bodies (Ida & Makina ll992t 
Ohtsuki & Emori 2000). Specifically, the Rayleigh distribu- 
tion 



2e 

P R (e\o- e ) = —r exp 



0-7, 



(32) 



is often considered, where cr e is the rms-value of the ec- 
centricity. Naively, one might expect that the distribution- 
averaged torque is just the torque expressions calculated listed 
in Table |2] averaged by Pit(eo\cr e o). However, this implies that 
the eccentricity distribution at a distance b is merely shifted 
in e, which is not the same as a gradient in cr e - the situation 
we consider here. 

Thus, we write cr e * cr e0 + f3bcr e Q, where [3 is the gradient 
with respect to cr e and b the dimensionless separation b = 
(a - ao)/flo- For t « 1 we can expand Equation d32i i with 
respect to b and write P R (e\cr e ) as P R (e\o- e0 ) + /3bPc(e\o-eo), 
where Pc is a correction term: 



Pc(e\o- e o) = 4e 



2 2 

e - °7o 



exp 



(33) 



The situation is illustrated in Fig. [3] At the reference point «o 
the velocity distribution is characterized by an rms-value cr e Q. 



For positive (3 the rms-value increases. The thin dashed line 
in Fig. |3] illustrates the case where cr e has increased by 10% 
to l.lCgrj, which corresponds to a distance b = 0.1 //3. Com- 
pared to the eccentricity distribution at b — 0, the distribution 
at b — 0.1 /yS has an excess of high eccentricity planetesimals 
and a deficit of low eccentricity planetesimals. To first order 
in b this change is represented by f}bPc{cr e \i). As a first order 
approximation, therefore, the distribution at b is well approx- 
imated by the superposition of PR(°~eo) and /3bPc(cr e o). 

The pbPc component is linear in b. The integration of this 
term over b thus behaves the same as the integration of the 
gradient in the surface density. To see this, recall that we lin- 
earized S = So(a/ao) a as S » Eo(l + ab) and found that the 
ab term gave rise to a torque cryv (see §[3]i. Analogously, 
integration of fibPc will yield a term fiy^Pc in the expres- 
sion for the dimensionless torque. The integration of the Pr 
component proceeds as before, except that /3 = since the 
eccentricity gradient is already included via Pc- Having ac- 
counted for the spatial distribution in this way, the resulting 
expression must yet be averaged over the Rayleigh velocity 
distribution at «o- Without loss of generality, we will consider 
a two-sided torque for which F(eo) - e Q 2 (see Tabled; the 
distribution-averaged torque then reads: Q 



<r 2 . s > 



(a £l ) 2 M, 



S 



deo 



= q p q d x (34) 

[%w + QTv] fgCfolfeo) +PjyP cCfolfeo) 
e 2 



The integration diverges for eo — > 0, which is because the 
shear-dominated regime is not covered in this work. There- 
fore the integration is cut off at the point where the Hill veloc- 
ity <?/, becomes 1 (see Equation Q). We find: 



<r 2s ) 



(a fio) 2 M p 



, 01 , % v + yv(a-2p) 

(2 log o- MJl )q p q d , (36) 



e0 



where cr^/, is the rms-eccentricity expressed in Hill units. 
Equation ( B6l l is very similar to the torque expressions ob- 
tained for the single-value power-laws. 

The above is merely a sketch for the inclusion of a velocity 
distribution. It is not complete, since a restriction on the incli- 
nation (i = e/2 or i <K e) is enforced. More realistically, the 
velocity distribution will be two dimensional and read, instead 
of Equation d32b : 



P R (i, e\cr e , cr,) 



Me 
ola 2 



exp 







2 




- 


\<r„ 







(37) 



Furthermore, when the close encounter torque is considered, a te rm — <5yv 
must be added to the terms within the square brackets Equation {34). This, to 
account for th e sc aleheight dependence. The <S-term will likewise propagate 
into Equation )36t . 

4 The integrals in Equation B41 evaluate to 



r 

■j t'i 

r 



e 2 



de z 



de 



r 0,4- 



2 exp 



y + 2 log ■ 



T 



(35) 



2(l +r + 21og^) 



r 2 

r e0 



where r(0, x) is here the incomplete gamma function, y at 0.577 the Euler- 
Mascheroni constant, and where we have expa nded with respect to ex, the 
lower cut-off of the integrations. In Equation j36\ we have only kept the 
logartihmic terms in the denominator and inserted e\ = CTgo/o-goj,. 
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In addition, we have in Equation (1361 1 neglected variations of 
the Coulomb factor, f A , with i and e. These may be come im- 
porta nt when accounting for a distribution average dlda et al.l 
1993). Finally, for a general treatment, an expression for the 
torque at arbitrary z'rj and eo is needed. A general expression 
is derived in Appendix IB. 3. 3 1 Therefore, the framework to 
(numerically) compute a truly Rayleigh-distributed average 
torque is in place. 

4. THE ROLE OF DIFFUSION IN CLOSE ENCOUNTERS 

Encounters with single planetesimals result in either inward 
or outward kicks dependent on the direction of the scattering. 
Migration therefore always has a stochastic (random) com- 
ponent. However, the inward and outward contributions are 
not equal. What we have calculated by the two-sided torques 
expressions of §0- the residual - is the systematic compo- 
nent. Which component will dominate depends on the ratio 
of the planetesimal mass over the planet mass, m/M p , and the 
lengthscale of interest. 

To quantify the migration rate due to stochastic motions, 
we calculate the diffusion coefficient, D[(Au^,) 2 ], which fol- 
lows from the change in th e diffusion rates of the para llel and 
perpendicular components (Binnev & Tremaine 2008): 



D[A Vi Avj] = ^ j£>[(A y ||) 2 ] - ^D[{Av ± ) 2 ]^ + yijD[{Av ± ) 2 ] 

(38) 



with 



D[(Av n ) 2 ]=47mv^—g A 

ir 

G 2 m 2 

D[(Av ± ) 2 ] = 4mv^-<J A - g A ) 



(39) 
(40) 



where g A = A 2 /(l + A 2 ). Thus, in our case we must calculate 



D[(Av<f>) z ] = 2nG 2 N m 2 { -^[3g A - f A ] + — [f A — g A ] \ . (41) 



We next perform similar operations as outlined in § 12.31 That 
is, we express n by the density function Pr, (Equation ( fl9l » 
and integrate over the semi-major axis. The steps are outlined 
in Appendix |C] We obtain 



D[(Av,) 2 ] 



3.5/ A 



l-5g A qdq P m 2 3 

flniin 



Mr, 



-H)"0- 



(42) 



Rather than the diffusion of (At^) we seek the diffusion in 
«q, D[(Aao) 2 ], which we refer to as the viscosity v. Since 
Aa = 2Aiu/Oo we have 



4D[(A^) 2 ] 14/ A - 6g t 



n 2 



a V m . (43) 



This is the diffusion coefficient for planets that results from 
the backreaction to the scattering of planetesimals. Contrary 
to the migration rate it does not depend on the exponents of 
a and ft, but it does involve the mass of the planetesimal (to). 
When to <s M p there are many encounters, whose individual 
kicks are small, resulting in a smooth migration rate. When 
to starts to approach M p , on the other hand, the importance of 
diffusive (random) motion increases. As a result, the migra- 
tion becomes increasingly 'noisy'. The critical lengthscale is 




Disk semi major axis, a 

Figure 4. Sketch of eccentricity profile in the self-regulated regime. A planet 
at ciq migrates in the direction of a less eccentric planetesimal disk (outwards 
in this sketch), interacting with planetesimal through close encounters over 
a width Aa. During its passage, the planet slightly excites the planetesimal 
disk from pre-encounter eccentricities e plc to post-encounter eccentricities 
£post. We solve for the jump in eccentricity Aeg = (e pos t ~ £pie)/2 and the 
corresponding eccentricity index /3 sr assuming that the relative increase in 
eccentricity is small, Aeo/eo <K 1. However, even though Aeo/eg « 1 the 
eccentricity exponent /3 S1 as Aeo/e^ can become large (s>l), affecting the 
migration rate. 



given by L* ~ JvT m i ai ~ Jm/M p yaQ. For L <sc L* diffu- 
fflc ' 



sive behavior will dominate. On scales L » L* the migration 
occurs smoothly. 



4.1. Comparison to Ohts uki & Tanakd ( 120031) 

The diffusion coefficient for planets (Equation ( l43l )) may be 
compared to the coefficient applica ble for an egual-ma s s plan- 
etesimal swarm as calculated by Ohtsuki & Tanakal d2003l) 
(their Equation [18]): 



VOT03 



: (N s a Q )(2 n m ) a Q ll, (44) 



ne h i h 



where A^ s is the column density, Iim = {qn/3) 1 ^ and 7rvs(j6) ~ 
0.3 for i = e/2. Expressing Equation ( l44l > in our notation, we 
find 

4/a m 



VOT03 



4' 9d9p M 1 



■a V k0 



(45) 



which is of the same magnitude as Equation (1431 1, Thus, 
a planet of mass M p interacting with a swarm of plan- 
etesimals of mass m diffuses at the same rate as the plan- 
etesimals do by interacting among themselves! Physically, 
the increase in cross section for encounters between planet 
and planetesimal due to its larger mass (M p ) is balanced 
by the decreasing kick (Aa) the planet receives. The latter 
scales as m/M p , while the cross section for encounters in the 
dispersion-dominated regime scales as cr oc b^ Q cc (GM p ) 2 
(e.g., Binnev & Tremainel 120081) . Thus, v oc (Aa) 2 cr stays 
constant. 

5. SELF-REGULATED PLANET MIGRATION 

The migration timescale (Equation d28l )) that we have de- 
rived assumes that the local distribution of eccentricity and 
surface density of the planetesimals is given by the power-law 
indices a and [} that characterize the protoplanetary disk on 
global scales. That is, any feedback of the planet on the disk 
that can change a and /3 locally is ignored. We will now relax 
this assumptions and consider the case where the protoplanet 
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slightly excites the motions of the planetesimals, increasing 
their eccentricity, and altering the local eccentricity profile. 

Specifically, we consider the configuration sketched in 
Fig. [4] where the eccentricity profile is steady in the frame 
of the migrating planet. The planet can migrate outwards 
or inwards (Fig. [4] depicts outwards migration). During its 
passage, which is defined as the time during which plan- 
etesimals interact through close encounters, planetesimals at 
a pre-stirring eccentricity e pie are excited to an eccentricity 
<? P ost- We assume that the eccentricity jump, Ae <K eo, where 
eo is the eccentricity at the reference radius. Moreover, we 
assume that the eccentricity changes gradually, which im- 
plies that a planetesimal experiences many (small) encoun- 
ters, before the planet has moved away from the interaction 
zone. Let the half- width of the interaction zone be denoted 
Aa * eoao <s ao- Then, we can expand Equation (HJ to obtain 
e ~ eo + ySeoAa/flo; therefore, Ae m f3e^ or 



Ae 



(46) 



This equation suggests that ft can become large, even though 
Ae /e <s 1. 

We seek to obtain an expression for the eccentricity jump, 
Ae. We may write 

A ^^S) eo ' (r - <<7 - } (47) 

where T* is the timescale to migrate locally over a distance 

migr ° J 

Aa = eo«o and T vs is the stirring timescale. Indeed, we should 
have that T\ <K T vs since we assume that Ae <k e . Phys- 
ically, this means that the planet migrates faster than it can 
excite the planetesimal disk. The local migration timescale 
T^j r is therefore just Equation d28l l multiplied by eo: 



12/ajS 



"o ' 



Id 1p 



(48) 



where we assumed that y t ot ~ Tci ~ _ 6/3/a is entirely due to 
the /3-dependence in y c \. The viscous stirring timescale in the 
dispersion-dominated regime is0 



2e 5 



On 



(50) 



With these expressions Equation d4Tb becomes 

Ae*-^. (51) 
24/3e q d 

Equating Equation ( Bit with Ae in Equation ( f46b we obtain 
the eccentricity power law index for self -regulated migration: 



c ll' 



1 



\-3/2 



(52) 



5 The viscous stir ring timescale is obtained from the relaxation timescale 
tChan drasekliar 1943; see also Ida & Mak ino1ll993) : 

T ys * T ch = , (49) 

Ann p (G N M p Y In A 

where n p is the density of perturbers. In our case n„ is the single protoplanet 
divided by its 'stirring volume' (2nao) X (2eoao) X (2ioOo). For v = egaoilo, 
'0 = «o/2, GnM+ = a^Q.^ and In A » /a/2, one obtains Equation i50i . 
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Figure 5. Migration timescales as function of eccentricity for a disk mass 
qj = 10~ 4 and a planet mass q„ = 10 -6 . Shown are the migration timescale 
in terms of the inverse orbital frequency corresponding to: (i) a fixe d di- 
mensionless torque value |y t otl = 10 (solid, black line; Equation (28)); (ii) 
the self-regulated migration scenario of §[5] where we s olve for the local ec- 
centricity gradient fl (dashed, bl ack line; Equation (53)); (iii) the one-sided 
torque only (gray line; Equation (30)); (iv) the type-I migration rate (dotted, 
horizontal line; Equation 13 11 ) using q CJ = 10 2 . O ur expression for T ST will 
break down below an eccentricity e" h (Equation (59)), indicated by a star. 

where we have switched to Hill eccentricities (see 
Equation ((5])). With Equation d52l ) we can solve for r and 
the (global) migration timescale in the self-regulated regime, 
T* r — T*. /eo: 



7/2 
f 



-.4 // 



(53) 



This much shorter timescale than Equation d28t can be at- 
tributed to the large yS sr value. For our fiducial param- 
eters j6 sl * 7. Nevertheless, the eccentricity jump Ae 
(Equation dBTl l) is only 3 x 10~ 3 , much less than eo. 

Figure [5] plots the migration timescale in the self-regulated 
limit, Equation (1531 1. as function of eccentricity for q p = 10" 6 
and q c i = 10 -4 . Also plotted is r m j gl of Equation d28l > 
for a constant value of y t ot = 10, the one-sided migration 
timescale (7Y S ), and the type-I migration timescale r t yp e -i 
(Equation dBTll). As can be seen from Fig.|5j self-regulated 
migration results in a migration timescale that is significantly 
shorter than Equation d28l ). especially at low and modest ec- 
centricities. In fact, it approaches the migration rate due to the 
one-sided torque; for e/, = 3-4 the self -regulated migration 
rate rivals that of as type-I. 

5.1. Preconditions for the self-regulated migration regime 

In the above analysis we have assumed that the planet exerts 
a modest feedback on the disk. For the self -regulated migra- 
tion mechanism to operate the planet mass can neither be too 
small (since then no feedback is present) nor too massive (too 
much feedback). If it is too small f3 sr (Equation ( l52b ) will be 
much lower than the global value of the eccentricity index /3. 
Thus, ^6 S1 | > \j3\ is a rough estimate of the precondition for the 
self -regulated mechanism to become feasible. 

Another assumption in the above analysis is that the jump 
in eccentricity gradient is smooth. More specifically, for the 
derivation of Equation d53b to be valid the timescale inequality 
relation 



Tsyn « K igl . « r Vi 



(54) 
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has to be obeyed. The first inequality ensures that the num- 
ber of scattering during the passage of the planet (TV ) is 
» 1 : the eccentricity of a planetesimal when it proceeds down- 
stream (in the frame of the planet) then increases gradually. 
The second inequality indicates that the total jump in eccen- 
tricity, Ae/<?o <k 1, is small: the planet is of a small enough 
size to only mildly excite the disk. 

The synodical timescale can be approximated as r syn = 
4n/3eoQ.{). Equation (l54l then reads 

„3 



An 

3e 



12./Aj3<7rf<2> 



2e 5 



(55) 



Converting to Hill eccentricities, eo = 
ranging, Equation (T55l > transforms into 



ei,{q p /3) 1 ^, and rear- 



6tt/ a « Q pd el l/2 « e 



where 



1/2 



(56) 



(57) 



is a combination of the planet and the disk masses (for the 
values of q^ and q p used in Fig. [5] 2pd = 0.25). With this 
notation, the second inequality of Equation d54l corresponds 
to 

e h » Q 2 pd , (58) 
whereas the first inequality becomes 
<6nf A 



» el 



2/11 



2.7 



2/11 



/ q P r 2/33 / \ 1/n 
\ 10 6 / \ 10 4 / 



(59) 

Both conditions must be satisfied. Note that when gpd ^ 1 
criterion Equation ( |58l vanishes since e/, is always larger than 
unity in the dispersion dominated regime. Physically, Qpd ex- 
presses the mass of the planet (the q p term) with respect to the 
planetesimal disk (the qd term). When Qpd is large, the planet 
has too much inertia to be affected by planetesimal scattering. 
However, even when Q pd <K 1 it is required that e/, > e* h to 
ensure a smooth power-law gradient in eccentricity - a pre- 
condition in the derivation of Equation ( T53l . 

6. DISCUSSION 

6.1. Migration regimes 

Using the results of § [5] we can identify three regimes, 
where the migration behavior is qualitatively different: 

1. The low planet mass regime, where the planet can be 
regarded as a test body that does not affect the structure 
(surface density, eccentricity) of the disk; 

2. A large planet mass regime, where the planet has too 
much inertia to experience sustained migration; 

3. An intermediate regime, where the planet self-regulates 
its migration. 

These regimes are indicated in Fig.|6]by the roman literals, I, 
II, and EH. The thick black lines denote the regime boundaries. 
The boundary between the first and third regime correspond 
to the criterion /3 sr = 1. The boundary between regime II 
and III follows from Equation (1581 . Contour lines give migra- 
tion timescales r m i gr in terms of the local orbital period. The 
migration timescales of planets in regime I are those given by 




Planet mass, q p 

Figure 6. Classification of migration regimes for a disk mass parameter qj of 
10~ 5 (top) and 10~ 4 (bottom) and contours (solid gray lines) of the migration 
timescales in terms of fl" 1 . For small planets (area I) the feedback of the 
planet on the disk is negligible and the planet migration timescales are given 
by Equation (28). Migration in regime II is suppressed due to the large mass 
of the planet. Intermediate-mass planets self-regul ate t heir migration (area 
III). The migration rate is then given by Equation )53t as long as e;, > e* h 
(Equation (35)). 

Equation ( 1281 ). whereas those in regime III obey Equation 031 
as long as <?/, > el (dashed line). In regime II, planetesimals 
are excited before migration becomes important. Note that 
the choices for the regime boundaries are somewhat flexible; 
in reality, the transition between the regimes will be broader 
than suggested by the sharp boundaries of Fig. [6] 

Both the migration rate as well as the position of the regime 
boundaries depend on the disk mass parameter qj (Equation 
([Toll). A large disk mass enlarges the fraction of the parameter 
space occupied by regime I; whereas a low disk mass enlarges 
that of regime EL The direction of the migration in regime I 
is defined by the exponents a and f3, characterizing the disk- 
wide power-law indices of the surface density and eccentricity 
profiles. In regime III, the migration direction is unspecified, 
as the migrating planet will self-adjust fi locally to ±/3 sr . That 
is, the direction depends on the history of the planet. 

For example, when a planet migrates outwards in regime I 
due to a negative gradient in the planetesimals' eccentricity, it 
may cross the boundary towards regime III, where it will con- 
tinue to migrate outwards at a faster pace. Similarly, a gradi- 
ent in qd (a local quantity) can also trigger a regime change. 
Generally, we find that migration timescales in regime I are 
long (unless qd is large) and that fast migration occurs in the 
self -regulated mode. However, regime I is important in setting 
the initial direction of the migration. 

The classification I, II, and III resembles the corresponding 
migration types for gas-driven migration. For type I, the feed- 
back of the planet on the disk is negligible. As shown in § 13.31 
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the expressions for the migration timescales show very similar 
scalings. Ty pe II migration is driven by diffusive motions of 
the gas (e.g., Gressel et al. 201 1) or solids. We showed in §[4] 
however, that for PDM these timescales are generally long. 
On the other hand, PDM in the type-Ill mode covers a large 
region of the parameter space (Fig.|6]l. 

6.2. Implications for N-body simulations 

Recently, state-of-the-art A^-body simulations investigating 
the migration behavior of a single planet embedded in a mas- 
sive p lanetesimal disk dKirsh et al.l 120091; ICapobian co et al.l 
201 1) exhibit an integration instability, i.e., a planet migrates 
o n timescales of the one- sided torque (Equation (1301 ). Since 
in iKirsh et all (|2009) and Capobia nco et al.l (1201 ll) the inter- 
actions operate mostly in the shear-dominated regime, scat- 
terings are strong and migration timescales shorfl Our work 
does not consider the shear-dominated regime; but similar to 
the above works, we discover a migration mechanism that 
is self-regulated. That is, the timescale to move away from 
the local stirring zone of the planet (7**., Equation (l48l >) is 
shorter than the viscous stirring timescale 7\, s . As long as the 
conditions remain in place (most notably, the disk in plan- 
etesimals must be mas s ive; se e Fig. [6j the migration does not 
stop. As in Kirs h et al.l ((2009), the migration mechanism does 
not specify a direction (inwards or outwards); this depends on 
some initial perturbation. Migration timescales down to 10 5 
times the local orbital period then seem quite viable. 

iV-body simulations are necessarily limited in terms of their 
dynamic range; due to their large numbers, it is often im- 
possible to resolve each planetesimal individually. Therefore, 
the superparticle concept is often employed, in which groups 
of planetesimals are represented by a single A^-body parti- 
cle, which effectively amounts to increasing their g ravitational 
mass m but keeping the surface d ensity constant ( Kirs h et all 
20091 iBromlev & K envon 2 0TPbl) . To save computational re- 
sources, one prefers a large amount of grouping (that is, a 
large m). But it is clear that a simulation involving too mas- 
sive superparticles will no longer accurately reflect the phys- 
ically system. This is immediately clear in the extreme limit 
when m approaches M p , in which case a planet is scattered 
by the debris! As we quantified in §2] a large m/M p ratio in- 
creases the importance of diffusive motions ('noise'), and we 
calculated the length scale over which diffusive noise can be 
expected to be dominant. 

Alternatively, A^-body simulations often account for dy- 
namical friction from a debris (small particle) compo nent an- 
alytically. In a recent study, Leinhardt et al.l d2009l) studied 
the behavior of an A^-body system interacting in coagulation 
and migration including a prescription for F<jf resulting from 
interaction with debris. However, in their work the debris is 
assumed to move on non-eccentric Keplerian orbits and the 
Keplerian shear is not accounted for - both effects render the 
dynamical friction force artificially large. (Indeed planets are 
seen to migrate very rapidly inwards!). The correct proce- 
dure should follow the lines of this work; that is, one must 
solve for the local velocity field. Clearly, we must generalize 
the calculations to include eccentric planets at arbitrary incli- 
nation and include encoun ters in the shear-dominated regime 
( IBromlev & Kenvonl201 lal) - issues that will be addressed in 
the future. The calculations presented in this work therefore 
have the potential to provide a significant boost in the accu- 

6 To see this one may substitute eo(ei, ~ 1) =s (q p ) 1 ^ in Equation j30\ . 



racy of A^-body simulations containing a debris component. 

7. CONCLUSIONS 

In this paper, we have employed detailed analytical and nu- 
merical calculations to obtain the net torque acting on a planet 
due to the recoil from gravitational interactions with planetes- 
imals in the dispersion-dominated regime. We have included 
both distant and close encounters and obtained the net migra- 
tion rate by summing the torques from the interior and the 
exterior disks. We list our conclusions: 

1 . While the magnitude of the migration rate is primarily 
determined by the local values of the surface density 
(So) and eccentricity (<?o), the direction is given by the 
local gradient in these quantities (a and /?) and by the 
Coulomb factor /a (Fig. |2j». Usually, the contribution 
from close encounters will determine the migration di- 
rection, unless /a (and by implication e) are low. 

2. The expressions for the migration timescale Equation 
(|28t due to planetesimal scattering display similarities 
to type-I migration (Equation OTTl). if one replaces the 
disk mass in planetesimals by that of the gas and the 
eccentricity by c s /Vk- Since the disk mass in solids 
is lower, the planetesimal-driven migration timescale is 
generally longer. 

3. Under certain conditions a much faster migration mode 
(rivaling that of type-I) is obtained when the feedback 
of the planet on the disk is accounted for. The planet 
then self -regulates the value of the eccentricity gradient 
to yS sr , which is a function of the local physical parame- 
ters (Equation d52l). Generally, |/J sr | » 1 and migration 
is quite rapid (Equation d28Ti). 

4. As function of the dimensionless disk mass qj (Equa- 
tion (flOli). planet mass q p = M p /M*, and planetesi- 
mal eccentricity e, we have identified three migration 
regimes (Fig. [6j representing: (I) low mass planets, 
for which disk excitation is negligible; (II) high-mass 
planets, too massive to migrate significantly; and (III) 
intermediate-mass planets, which exert a mild feedback 
on the disk and migrate in the self -regulated mode. 
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APPENDIX 

A. TORQUE DENSITY AND MIGRATION RATE FOR DISTANT ENCOUNTERS 

A planet in a gas or particular disks exerts a torque on the disk at discrete positions (resonances). These resonances are 
identified by integer numbers I, m, which refer to the corresponding Fourier modes. For given m, the strongest resonance occurs 
where I = m, for which the resonance condition (a relation between m and the distance to the planet b) reads: 



+ -) = £2o, 
\ ml 



(Al) 



where e = ±1 = sgn(b) refers to an inner or outer Lindblad resonance and s — refers to a co-rotation resonance. 

Here we consider Lindblad resonances. For m » 1 the resonances start to overlap and a continuum tre atment is po ssible. 
Therefore, a torque density, dT/dr, can be defined (Goldr eich & TremaineHl980h . We follow the notation of (Ward 1997), who 
defines: 



dY 

dr 



K [r] 



(A2) 



where ip m is an expression that inv olves the computation of Laplace coefficients, and k is the epicyclic fre quency at r. It is 
important that the RHS of Equation ( IA2b is evaluated at resonance, whose condition is giv en by Equation ( IA1 b . 
We redefine b such that it becomes dimensionless, b — > b/ao <K 1, and expand Equation (lAlb with respect to b: 



— -i.fi-*). 

1 31*1 V, 4/ 



(1 + b) 3 ' 2 ■ 

where e = sgn(£>). IWardl d 1997b has obtained the function ip 2 n to first order in m (or b). He finds 



i[, 2 (b) = [2tf<>(2/3) + Kd2/3)Y 



| b ( *o(2/3) + 5^(2/3) 1 
2 12*0(2/3) + *! (2/3) J 



6.35(1 + 1.255/7) 



(A3) 



(A4) 



Thus, the four ^-dependent terms in Equation ( IA2b expand as follows: ii 2 (b) as in Equation JA4t : l/r » (1 + b)/ao; m 4 * 
(2/3/>) 4 (l - b); and the epicycle frequency as (Qo/k) 2 a 1 + 3b since K(r) = Q.(r) in a Keplerian disk. Collecting these terms gives 
the torque density to first order in b: 



dr 2VZ fl 4o2 

db * Sgn(b) W [2Ko(2/3) + ^i( 2 / 3 )] d + 2.255fe), 



(A5) 



where K v (x) is the modified Bessel function of second kind of order v (Note again that in this equation b is nondimensional). To 
zeroth order (£ = So) therefore: 



= sgn(/7)(M pfl 2 Q 2 )^|^ [2*o(2/3) + *i(2/3)] 2 * 2.51(M p a^) sgn(fc)^, 



(A6) 



where we substituted q p = Mp/M* and qj = a^Ho/M*. 

The torque from the disk on the planet has the opposite sign as Equation (IA5t . When we consider one only side of the disk, 
the leading term will not vanish. For example, the torque from the outer disk on the planet is: 



ls-di 



M^flo^o) 2 



db 



dTo 
db 



2 5 q p q d 
35 e 3 



+0.836M, 



4 



(A7) 
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(where the sign is positive for an interior disk). This expression is consistent with the migration rate der ived in Equation (0. 

To zeroth order, the integration over both sides of the disk vanishes due to symmetry of Equation (IA6t with respect to b. In the 
next order expansion of Equation (IA5t . however, nonzero contributions originate due to: 



1. The curvature term in Equation ( 1A51 >: and the gradient in surface density a. These terms render the integrand odd (ocl/fo 3 ) 
and together contribute a factor 



(2.26 + a)b = -3e (2.26 + a) |r Is _ di | (A8) 



2. The gradient in eccentricity. This causes the transition between close and distant encounters to shift by a value b = fieo 
(see Fig.Q]): from -eo < b < eo (when f3 = 0) to -eo + /6e ( 2 < eo + jSe^. The corresponding net torque due to this shift is the 
zeroth order torque density at b — eo times twice the width of this shift 



db\— j xlfietx^—] [b = eo] = 6/teoiri s -dil (A9) 



Summing the two terms gives the leading term of the torque on the planet when accounting for both sides of the disk: 

r 2s -di = -3e (2.26 + a - 2/3W ls -J = M p (ao^o) 2 q d q P ~ 5 66 " " — ■ (A10) 

e o 

B. CALCULATION OF INTEGRALS FOR CLOSE ENCOUNTERS 

Let us write the Equation (l2"TT i in nondimensional form: 



a 2 M p Q. 2 (floQo) 4 

= q d q p X / c i(a,j8) 



J™ 3 "" da , - (a 
— a 2 P R A- 



1+a v^V 2 

U <l> y K0 



(Bl) 



where we used (G^M*) 2 = (fl^Qg) 2 . The integral in the square brackets is defined I c \(a,/3) and must be computed. 

The integration is over the range in semi-axes a where planetesimals are able to cross the planet's orbit ao. In the above v - |v| is 
the velocity of an unperturbed body at a — ao relative to the circularly-moving planet and the azimuthal component of v. Both 
v and Pr : , the probability density of planetesimals at ao (see below), are functions of the semi-major axis a of the planetesimals. 
In the following sections we will obtain expressi ons for v and Pr z , respectively. For aesthetic purposes we will express these 
quantities, as well as the other terms in Equation (IB lb . as function of 6 - the true anomaly of the Kepler orb it at the point where 
it intersects the planet - rather than a. A perturbation analysis in eo then allows us to compute Equation ( IBU analytically. 

B.l. Expressions resulting from the Kepler orbit 
A body traveling on a Kepler orbit obeys the relation 

r- a(1 ~ e2 \ (B2) 
1 + ecosd 

where r is the radial coordinate in the orbital plane of the particle, a the semi-major axis, e the eccentricity, and 9 the true anomaly 
specifying the instantaneous position of the particle. This orbital plane is inclined by an angle i with respect to the equatorial 
plane - the plane in which the planet moves. The particle intersects the planet's circular orbit at r = ao, or at a true anomaly 8 
that obeys the relation: 

1 + ecosf? 

a = a 9 = a — =— . (B3) 

(1 -e z ) 

Thus, there is a one-to-one relation between the semi-majo r axi s a from which the planetesimal originates and the true anomaly 
9 at which it crosses the planet at ao- We will use Equation (IB3b to switch the integration variable to 9: 

da aoesin9 

Te = ~T^- (B4) 



Next, we express the Keplerian velocities also as function of 9. Kepler's second law gives the azimuthal velocity, Vg = l/r = 
yja(l - e 2 )p/r with / the angular momentum and p = Gn(M+ + M p ) « G^M*. The radial velocity can be obtained from energy 
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Planet, e = 

a = ao/(l + <?) 
a = ao 

a = flo/(l — e ) 



Figure 7. Sketch of sample trajectories in the orbital plane. A planet moves on a circular orbit of radius ciq (blue circle) and interacts with planetesimals that are 
characterized by semi-major axis a and eccentricity e = 0.3. For fixed e the range in a that crosses the planet's orbit is ao/(i + e) < a < ag/(l - e). The periapsis 
(8 = 0) of the orbits are, for clarity, situated on the positive X axis. The planet interacts with the outermost planetesimal swarm (gray, solid circle) at their 
periapsis (8 = 0) and with the innermost swarm (gray, dotted circle) at their apoapsis (6 = n). For planetesimals of intermediate a, e.g., those of a = ag (black 
circle), 8 is given by Equation )B3t . At the interacti on p oint ( see inset) the planet moves along the 8 direction at the Keplerian velocity 14o- The components of 
the planetesimal velocity V are given by Equations )B7t and )B8t . For the calculation of the dynamical friction force it is the relative velocity v = V - Vjo that 
matters (red arrow). 



conservation. We evaluate these expressions at the interaction point, i.e., at r — ao and a/ao given by Equation ( IB3b : Q 



G N M+(l-e 2 ) a r 

= Vfco VI + e cos 0; 

a ao 



V r = V a - 



«0 

esin6> 



Vl + ecos0 



(B7) 
(B8) 



where V^o = -^u/ao is the orbital velocity at the interaction point (the Kepler velocity at which the planet moves). In disk (i.e., 
cylindrical) coordinates (R, (p, z) the velocities at the interaction point (R, z) = (ao, 0) at arbitrary tp become: 



V R = V r ; 

V = V e cos z; 

V z = Ve sin 

and the relative velocity vector, written in terms of 0, reads: 



v = V-V„ 



( V R 



= V k o 



e sin 9 



vr 



+ e cos v 
cos i + e cos 8-1 



sin i Vl + e cos ( 



(B9) 
(BIO) 
(Bll) 



(B12) 



7 The expression for the radial velocity expression is perhaps difficult to and Equation m%\ is retrieved, 

see at first glance. Energy conservation gives — 



}_ (V 2 „2, __ M_ •E. = JL(i-^.\= V *0 / t+Zecosg + e 2 
2^ r 9 2a a a Q \ la) 2 \ l + ecos0 

where we used Equation )B3t . Then, inserting Equation )B7t for Vg: 



(B5) 



l+2ecose + e 2 (l + ecos9) 2 \ 

— (Bo) 

1 + e cos 8 1 + e cos 8 J 
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B.2. The projection operator Pr z 

The orbit of a planetesimal is additio nally determin ed by o>, the angle of periapsis. Transforming from (r, ff) to Cartesian 
coordinates in the equatorial plane gives (Beutler2005): 



x = rg cos(a> + ff); 
y-rg sin(w + 9) cos 
z-r g sin(w + ff) sin i\ 



(B13) 
(B14) 
(B15) 



where we have assumed, without loss of generality, that the line of nodes is directed along the x-axis. From these equations we 
obtain the projected radius on the equatorial plane R: 



R 



= V'" 2 _ z 2 = r ^1 - sin 2 (w + ff) sin 2 i. 



(B16) 



For the problem considered here, R and z are the principal variables. The density function Pr z is defined such that P Rz dRdz gives 
the probability of finding the particle within [R, R + dR; z, z + dz]. To obtain Pr z , we assume that the phase angles t (mean anomaly) 
and a> are randomly distributed, P tM = Q. a /(2n) 2 , where Q fl denotes the orbital frequency corresponding to semi-major axis a. 
Converting variables then gives: 

d(t, o>) 



Prz 



d(R,z) 



PtA 



(B17) 



where 8(t, o>)/8(R, z) is the Jacobian of the transformation. 

We may proceed to use Kepler's equation to relate the mean anomaly t to 9. Here, howev er, we consider a specific case 
where P Rz is only evaluated at (R,z) = (ao,Q). Let this density be denoted P Rz . From Equation dB15l ) it can be seen that z = 
corresponds to either a> — -9 or co — n - 9. Therefore, Pr z is a function of one variable only, say 9. Formally, we define 



Prz = Prz(.z = 0) 



d(t, oS) 



d(R,z) 



P tt0 [6(0} = -ff) + 6(0) = 7T - ff)] ., 



(B18) 



where 6(x) is the Dirac delta function. Using Equation dB 16b and Equations ( IB13t - dBT5b and anticipating that the matrix elements 
of the Jacobian will be evaluated at rg = qq and o) — -9 or o> — n - ff. 



3z_ 
do) 



= r g cos(w + 9) sin i — > aq sin i 



dz dr 

— - rg cos(<y + 9) sin i H sin(w + 9) sin ; 

89 89 







8R _ z(8z/do)) 

80) y/ r 2 + z 2 

8R _ r(8r/89) - z(8z/89) 
89 ~ 



Vr 2 - z 2 



8r_ 

89 



8R 



flo sin i 



= V r 



(B19) 
(B20) 
(B21) 

(B22) 



where '— >' indicates we have evaluated the matrix elements at z — (or o) — -9) and R = ao- The Jacobian, evaluated at 
(R, z) = (ao, 0) then reduces to 

rift c,A rtv 1 

(B23) 



(B24) 



8(t, o)) 




8R dz 


1 


8(R,z) 


R=ao;z=0 


8t do) 


R=a -z=o V r d sin i 



an expression that only involves 9 as a variable. Substituting, we obtain for the density function Equation ( IB 18b 

Cl a (9) 



Prz(0) 



2n 2 V r ao sin i 



Substituting Equation (|B8t for V r and Q.„ = Qo( fl fl/ fl o) 3 ^ 2 we nave 



alP Rz (9) 



a„ 



Vl + ecos# 



2n 2 e sin 9 sin i \ao 



■3/2 



(B25) 



Equation ( IB24l i can be validated numerically, for example by distributing particles on a Kepler orbit characterized by fixed orbital 
elements a, e, and i but random o) and t. The midplane number density corresponding to the reference radius ao, «mid(ao[#]), is 
then «mid(ao) = ^Prz/^^o- The additional factor of 2 takes care of the fact that there are two ^-solutions for a given ao- 



B.3. Integral evaluations 

Using Equation (IB4b to transform coordinates to 9 and Equation (IB25b for c&Prz, I c \(a,/3) (see Equation (IB It ) can be expressed 
solely as a function of 9: 



I a = In ft 



I 



d0 



Vl +ecos 9 la g \ a 1/2 v^Vj 



2n 2 sin z'(l - e 2 ) \aoj 



r- 



(B26) 
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where ag, v$ and v are given by Equa tions ( IB3t and (IB 121 >. Furthermore, e and i are also functions of 9 by virtue of the gra- 
dient (Equation (Q]i). Since Equation ( IB26b cannot be solved algebraically, we expand it in e to obtain a closed-form solution. 
Concerning the inclination, we will consider two cases: (i) i = e/2 <K 1; and (ii) i « c « 1. 

B.3.1. The equilibrium solution, i = e/2 
Substituting i — e/2 in Equation ( IB26I ) and approximate <?, assuming e <K 1: 

e(ff) w e + y6<?o cos 0. (B27) 

Next, we expand Equation ( lB26b around eo, accounting only for terms of <9(e Q 3 ) and 0{e^ 2 ). The integrand of Equation ( 1B261 
then becomes: 

( 44cos(0) - 12cos(36») _ 3 22a - 66/? + 2(8a - 24/? + 1) cos(2<9) - 3(2a - 6/3 + 3) cos(4<9) - 9 _ 2 | 

" + ■ ■ -775 e n > . (B28) 



(2^(cos 2 (0) + 4sin 2 (0) + l) 2tt(cos 2 (6») + 4sin 2 (6») + l) j 

The e Q 3 term is symmetric and evaluates to when the full range of 9 is considered. When we integrate only over one side of the 
disk (e.g., -7r/2 < < n/2) this term determines the one-sided torque: 

ls - d V5,re 3 el 

The e 2 term in Equation ( IB28b does not vanish after integration over < 9 < 2n and evaluates to 

4 V2 (g (-1) (12. - 36/? + 5) + K (-§) (-12a + 3^- 11)) -0,66+ 1.97(^-3/?) , 
/ 2s _ c] = — ; A* A- ( B3 °) 

where £(x) and K(x) are complete elliptic integrals of the first and second kind: 

c/0 r 71 ^ 2 / 

: £(*) = aWVl -xsin 2 6». (B31) 
o Vl - Jtsin 2 <9 -'o 

B.3.2. The thin disk case, i <k e 

Next we consider a case of a thin disk (i <k e). A thin disk is applicable when the planet interacts with the planetesimals in the 
shear-dominated regime. Note that we have assumed in the main text that the dispe rsion-dominated regime a pplies, for which 
i/e « 0.5 is expected, but a situation where i <K e can still emerge as a transient state (Rafikov & Slepian 2010). 

We consider the situation where the inclination also obeys a power-law: 

/ a \ S 

i(a) — ;'o — ~ io + <5i'oeo cos 9. (B32) 
\ao) 

We follow the same procedure as above, expanding / c i first in terms of z'o and then in terms of e . Subsequently, integration over 
9 gives: 

M<e) _ 4 /a _ L27 /a , R „> 



12 [-2K(-3) -K($) + 2E(-3) + 4E (|)] (a -2/3 -5)- 14 [2K(-3) + K (\)\ + 44£ (f) + 22£(-3) 

9jre jo A 

1 63 + 2.2RI7* - 2/? - #) 

-/a. (B34) 



Wo 

for the one- and two-sided integrals, respectively. 

B.3.3. General case 

Within the limits of our assumptions (inclin ation s and eccentrici ties la rger than the Hill eccentricity, etc) we consider an even 
more general case. After inserting Equation dB27| ) and Equation dB32| ) for e{9) and i(9), respectively, we define ig = (eo and 
expand I c \ in terms of eo. The resulting expression (a function of f and 8) can be integrated. This procedure gives: 

/u = : _j£ (B35) 

7 2s = V ^ , * ' V /a (B36) 

97re 2 V4<T 2 + 1 (4<T 5 + 5<T 3 + £) 
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where 

Ai=4^(12a- 12/3- 246- l) + £ 2 (60or- 36/3- 1445 + 7)+ 12a -24/?- 125+11 (B37) 
A 2 = 8^ 2 (-3« + 3/3 + 66 - 2) - 6a + 12/3 + 66 - 7 (B38) 

It can be verified that these formula recover the expressions for the equilibrium regime (where £ = 1/2 and 6 = /3) and the z <k e 
regime (where f « 1), respectively. 

C. CALCULATION FOR THE DIFFUSION INTEGRALS (CLOSE ENCOUNTERS) 

Equation ( l4TT i gives the rate of change in Ai>^: 

f nv \ n ) 

D[(Au^) 2 ] = 2nG 2 N m 2 1 -f[3g A - / A ] + -[/a - 2a] ■ (CI) 

The number density «, and v are functions of disk radi us a and Equation dCU must accordingly be converted in an integration 
over a. Following a similar procedure as described in § 12.31 we write: 

D[(Avt) 2 ] = 2nG 2 N m J da P R X(a) {^j j(3g A - / A )^ + • (C2) 

After inserting Equation (Q} for X(a) we split the integral, D[(A^) 2 ] = Do(h + h), with Do the dimensional part given by 

and /i and 7 2 dimensionless integrals given by 

, l+a ,,'2 



D = -4 — - = a 2 0; (C3) 



da'P'J-j ± (C4) 
h=2n{f h -g K ) J da'P' Rz ^j ^ (C5) 

where primes denote normalized quantities. 

The procedure to evaluate integrals Equation ( lC4t is the same as in Appendix[B] First, we chang e variables to 9 via Equations 
( lB3t and dB4l . Then we insert Pr z (9) and the velocity fie ld, \( 9), as given by Equations ( IB25I ) and ( IB12l i. respectively. We 
further assume equilibrium, i = e/2, insert e{9) (Equation (IB27I )) and expand I\,h in eo, keeping only the highest order term. 
Subsequently, we derive: 

8V2( /A -3 gA )(g(-i)-g(-i)) 0.98(3 gA -/ A ) 

n = 71 ~ 5 

3e n e 

8V2(/A- gA )^(-|) 4.4(/ A - gA ) 
'2 = 5 ~ 5 (<~ /) 

h + h« 3 - 5fA - 2 h5gA . (C8) 
p 

o 

Note that there is no dependence on a or ft as the highest-order terms do not cancel. 



